
#include <fstream>
#include <queue>
#include "mesh.h"

/// A functor for computing plane equations
struct PlaneEquation {
   template <class Facet>
   typename Facet::Plane_3 operator()(Facet& f) {
      typename Facet::Halfedge_handle h = f.halfedge();
      typedef typename Facet::Plane_3  Plane;
      return Plane(h->vertex()->point(),
                   h->next()->vertex()->point(),
                   h->next()->next()->vertex()->point());
   }
};

/// A functor for computing vertex normals 
struct VertexNormal {
   template <class Vertex>
   void operator()(Vertex& v) {
      Mesh::Halfedge_around_vertex_const_circulator 
            hbegin = v.vertex_begin(), h = hbegin;
      int n = 0;
      Vector3 sum(0,0,0);
      do {
         if (not h->is_border()) {
            sum = sum + h->facet()->plane().orthogonal_vector();
            n++;
         }
      } while (++h != hbegin);
      if (n>0 and sum*sum>0) {
         sum = sum/n;
         sum = sum/sqrt(sum*sum);
      }
      v.normal(sum);
   }
};

struct compareuv {
   bool operator() (const Vertex_handle &u, const Vertex_handle &v) const {
      return u->d > v->d;
   }
};

Mesh::Mesh() {
   
   _nvertices=_nfaces=0;
   _radius = 1.0;
}

Mesh::Mesh(const char* filename) {
   
   _nvertices=_nfaces=0;
   _radius = 1.0;
   loadShape(filename);
}

Mesh::~Mesh() {
   
}

void Mesh::loadShape(const char *filename) {
   
   cout<<"Loading... ";
   
   std::ifstream fin(filename);
   fin>>*this;
   fin.close();
   auxiliarData();
}

void Mesh::saveShapeOFF(const char *filename) {
   
   if (_nvertices==0 or _nfaces==0)   return;
   
   std::ofstream fout(filename);
   fout<<*this;
   fout.close();
}

Mesh* Mesh::load_off(const char *filename) {
   
   std::fstream off_file(filename, std::ios_base::in);
   if (off_file) {
      Mesh *p = new Mesh(filename);
      return p;
   }
   return 0; 
}

void Mesh::auxiliarData() {
   
   if (not this->is_pure_triangle()) {
      cout<<"### *** Attention . Attention . Attention . Attention *** ###"<<endl;
      cout<<"This example requires the mesh model to have triangle faces only!!!!."<<endl;
      return;
   }
   
   this->normalize_border();
   _nvertices = this->size_of_vertices();
   _nfaces    = this->size_of_facets();
   
   cout<<endl<<"Loading Done!"<<endl;
   cout<<"Vertices Number: "<<_nvertices<<endl;
   cout<<"Triangles Number: "<<_nfaces<<endl;
   if (this->size_of_border_edges() > 0)
      cout<<"The mesh has border edges "<<endl;
   else
      cout<<"The mesh hasnt border edges "<<endl;
   cout<<"KBytes used by the model: "<<this->bytes()/1000<<endl;
   cout<<"KBytes reserved by the model: "<<this->bytes_reserved()/1000<<endl;
   
   /// Compute the normal vector of each face
   std::transform(this->facets_begin(), this->facets_end(), 
                  this->planes_begin(), PlaneEquation());
   
   /// Compute the normal vector of each vertex
   std::for_each(this->vertices_begin(), this->vertices_end(), VertexNormal());
   
   /// Build the array containing the index of faces
   buildFacesIndices();
   
   /// Computes the limits of the model
   limits();
   
   /// Centers and scales the model
   //~ centerAndScale();
}

void Mesh::buildFacesIndices() {
   
   cout<<"buildArrays begin ... "<<endl;
   faces.resize(3*this->size_of_facets());
   // Build the array of faces
   unsigned iif=0;
   for (Facet_iterator fi=this->facets_begin(); fi!=this->facets_end(); ++fi, iif+=3) {
      Halfedge_around_facet_circulator j = fi->facet_begin();
      
      const unsigned &id0 = std::distance(this->vertices_begin(), j->vertex());   ++j;
      const unsigned &id1 = std::distance(this->vertices_begin(), j->vertex());   ++j;
      const unsigned &id2 = std::distance(this->vertices_begin(), j->vertex());
      faces[iif  ]=id0;
      faces[iif+1]=id1;
      faces[iif+2]=id2;
   }
   cout<<"... done! "<<endl;
   
   cout<<"... done!"<<endl;
   // Insert id at the halfedges
   unsigned int hii=0;
   for (Halfedge_iterator hi=this->halfedges_begin(); hi!=this->halfedges_end(); ++hi, ++hii) {
      hi->hedgeid = hii;
   }
   cout<<"... done!"<<endl;
   
   cout<<" begin ... "<<endl;
   // vertex ids
   unsigned iiv=0;
   for (Vertex_iterator vi = this->vertices_begin(); vi!=this->vertices_end(); ++vi, ++iiv) {
      vi->id(iiv);
   }
   cout<<"... done! "<<endl;
}

void Mesh::limits() {
   
   Vertex_iterator vit = this->vertices_begin();
   const Point3 &p = vit->point();
   for (int i = 0; i < 3; ++i) { mincoo[i] = maxcoo [i] = p[i]; }
   
   while (vit != this->vertices_end()) {
      const Point3 &p = vit->point();
      for (int i = 0; i < 3; ++i) {
         if (p[i] < mincoo[i])      mincoo[i] = p[i];
         else if (p[i] > maxcoo[i]) maxcoo[i] = p[i];
      }
      ++vit;
   }
   
   /// building the minimum bounding cube 
   Point3 WorldpMin(mincoo[0],mincoo[1],mincoo[2]);
   Point3 WorldpMax(maxcoo[0],maxcoo[1],maxcoo[2]);
   
   cout<<"pmin: "<<mincoo[0]<<" "<<mincoo[1]<<" "<<mincoo[2]<<endl;
   cout<<"pmax: "<<maxcoo[0]<<" "<<maxcoo[1]<<" "<<maxcoo[2]<<endl;
   
   Vector3 size(WorldpMin, WorldpMax);
   size = size*0.5;
   _center = WorldpMin + size;
   _radius = sqrt (size * size); // vector lenght
   size = Vector3(_radius, _radius, _radius);
   WorldpMin = _center - size;
   WorldpMax = _center + size;
   
   mincoo[0]=WorldpMin[0]; mincoo[1]=WorldpMin[1]; mincoo[2]=WorldpMin[2];
   maxcoo[0]=WorldpMax[0]; maxcoo[1]=WorldpMax[1]; maxcoo[2]=WorldpMax[2];
}

void Mesh::centerAndScale(Number ScaleFactor) {
   
   assert(ScaleFactor > 0);
   
   Number dimv[]={maxcoo[0]-mincoo[0], maxcoo[1]-mincoo[1], maxcoo[2]-mincoo[2]};
   assert(dimv[0]>0 and dimv[1]>0 and dimv[2]>0);
   
   // Determine the scale factor
   Number Scale;
   if (dimv[0] >= dimv[1] && dimv[0] >= dimv[2]) Scale = ScaleFactor/dimv[0];
   else if (dimv[1] >= dimv[0] && dimv[1] >= dimv[2]) Scale = ScaleFactor/dimv[1];
   else Scale = ScaleFactor/dimv[2];
   
   // Determine the translation vector
   Number transv[]={mincoo[0]+maxcoo[0], mincoo[1]+maxcoo[1], mincoo[2]+maxcoo[2]};
   transv[0] *= 0.5f;
   transv[1] *= 0.5f;
   transv[2] *= 0.5f;
   
   Vertex_iterator vit = this->vertices_begin();
   while (vit != this->vertices_end()) {
      Point3 &p = vit->point();
      p = Point3(p[0]-transv[0], p[1]-transv[1], p[2]-transv[2]);
      p = Point3(p[0]*Scale, p[1]*Scale, p[2]*Scale);
      ++vit;
   }
}

Point3 Mesh::getMidPointBetweenTwoFaces(uint idSelectedFront, uint idSelectedBack) {
   
   Facet_iterator ffi = facets_begin();
   advance(ffi, idSelectedFront);
   Mesh::Halfedge_around_facet_circulator fh = ffi->facet_begin();
   const Point3 &p0 = fh->vertex()->point();
   const Point3 &p1 = fh->next()->vertex()->point();
   const Point3 &p2 = fh->next()->next()->vertex()->point();
   Facet_iterator bfi = facets_begin();
   advance(bfi, idSelectedBack);
   Mesh::Halfedge_around_facet_circulator bh = bfi->facet_begin();
   const Point3 &q0 = bh->vertex()->point();
   const Point3 &q1 = bh->next()->vertex()->point();
   const Point3 &q2 = bh->next()->next()->vertex()->point();
   return midpoint(centroid(p0,p1,p2), centroid(q0,q1,q2));
}

void Mesh::computeLaplacianParameters(Number lfalpha, Number largest_dvj) {
   
   /// Compute the Laplacian parameters
   cout<<"Compute Laplacian parameters ... "<<endl;
   Vertex_iterator vi = this->vertices_begin();
   while (vi != this->vertices_end()) {
      if (vi->dClosestJoint > largest_dvj) largest_dvj = vi->dClosestJoint;
      ++vi;
   }
   
   cout<<" >>>largest_dvj... : "<<largest_dvj<<endl;
   
   Number lfscalek = lfalpha*largest_dvj*largest_dvj;
   vi = this->vertices_begin();
   while (vi != this->vertices_end()) {
      setVoronoiArea(vi);
      Number d2 = vi->dClosestJoint;
      assert(d2*d2 != 0);
      Number b = (lfscalek*vi->sapb)/(d2*d2);
      if ((b + vi->sapb) == 0) {
         cout<<"b + vi->sapb: "<<b<<" "<<vi->sapb<<endl;
      }
      assert((b + vi->sapb) != 0);
      vi->lfdenom = 1.0/(b + vi->sapb); // Laplacian filter denominator
      vi->lfkterm = b*vi->lfdenom; // Laplacian filter constant term
      ++vi;
   }
   
   ofstream saida("saida.txt");
   vi = this->vertices_begin();
   while (vi != this->vertices_end()) {
      saida<<vi->lfdenom<<" : "<<vi->lfkterm<<" >>> "<<vi->sapb<<endl;
      vi++;
   }
   cout<<"Compute Laplacian parameters ... end"<<endl;
}

void Mesh::setVoronoiArea(Vertex_handle v) {
   
   ofstream saida("saida.txt");
   v->VA = 0.0;  v->sapb = 0.0; // Initialize sum variables
   Halfedge_handle hevc = v->halfedge(), hhbegin = hevc;
   do {
      // Compute the alpha+beta for edge hevc
      setAlphaPlusBeta(hevc);
      
      // Compute the summation term for edge hevc
      Vertex_handle vi = hevc->vertex();
      Vertex_handle vj = hevc->opposite()->vertex();
      const Point3 &_pi = hevc->vertex()->point();
      const Point3 &_pj = hevc->opposite()->vertex()->point();
      const Point3 &_pl = hevc->next()->vertex()->point();
      Point3 pi(_pi[0],_pi[1],_pi[2]);
      Point3 pj(_pj[0],_pj[1],_pj[2]);
      Point3 pl(_pl[0],_pl[1],_pl[2]);
      Triangle3 T(pl,pi,pj);
      
      Vector3 pipj = pj-pi;
      double dot = pipj*pipj;
      
      // Accumulate the Voronoi area
      assert(dot != 0.0);
      if (CGAL::angle(pl,pi,pj) == CGAL::OBTUSE)
         v->VA += sqrt(T.squared_area())/2.0;
      else if (CGAL::angle(pi, pl, pj) == CGAL::OBTUSE or CGAL::angle(pi, pj, pl) == CGAL::OBTUSE)
         v->VA += (sqrt(T.squared_area()))/4.0;
      else
         v->VA += (hevc->getAlphaPlusBeta())*(pipj*pipj)/8.0;
      
      
      //~ if (v->id() == 1978) {
         //~ cout<<"     ::"<<v->sapb<<" : "<<v->VA<<" dot: "<<dot<<endl;
      //~ }
      
      // Sum cot values of opposited angles to current halfedge
      v->sapb += hevc->getAlphaPlusBeta();
      
      // Store information of incident halfedges of v
      v->neighbors.push_back(vj);
      v->cots.push_back(hevc->getAlphaPlusBeta());
      
      // Next halfede
      hevc = hevc->next()->opposite();
   } while(hevc != hhbegin);
   
   if(v->sapb == 0 or isnan(v->sapb)) {
      cout<<"id: "<<v->id()<<" point: "<<v->point()<<endl;
      Halfedge_handle hevc = v->halfedge(), hhbegin = hevc;
      do {
         cout<<"    "<<hevc->opposite()->vertex()->point()<<endl;
         cout<<"    "<<hevc->opposite()->vertex()->point() - hevc->vertex()->point()<<endl;
         hevc = hevc->next()->opposite();
      } while(hevc != hhbegin);
      
      assert(false);
   }
}

// Calcula la cotangente del angulo formado por los vectores A y B
inline Number cot(const Vector3& A, const Vector3& B) {
   
   if ((A*A)*(B*B) - (A*B)*(A*B) <= 0.0) {
      return 0;
      //~ cout<<" A, B:"<<A<<" "<<B<<endl;
      //~ assert(false);
   }
   
   if(isinf((A*B)/(std::sqrt(abs((A*A)*(B*B) - (A*B)*(A*B)))))) {
      cout<<" A, B:"<<A<<" "<<B<<endl;
      assert(false);
   }
   
   //~ assert((A*A)*(B*B) - (A*B)*(A*B) >= 0.0);
   Number AA = A*A;
   Number BB = B*B;
   Number AB = A*B;
   return AB/sqrt(abs(AA*BB-AB*AB));
   //~ return (A*B)/(std::sqrt(abs((A*A)*(B*B) - (A*B)*(A*B))));
}

void Mesh::setAlphaPlusBeta(Halfedge_handle he) {
   
   //~ cout<<" Halfedge_handle:    "<<
   
   Number alpha, beta; double epsilon = 1.0e-16;
   const Point3 &pi = he->vertex()->point(); // the vertex vi
   const Point3 &pj = he->opposite()->vertex()->point();  // the vertex vj : opposite to vi
   const Point3 &pb = he->next()->vertex()->point();      // the vertex vi-1
   const Point3 &pn = he->opposite()->next()->vertex()->point(); // the vertex vi+1
   
   // cuidando das bordas
   if (he->opposite()->vertex()->id() == he->next()->next()->vertex()->id())
      alpha = cot(pi-pb, pj-pb);
   else alpha = 0.0;
   if (alpha < epsilon) alpha = 0.0;
   
   if (he->vertex()->id() == he->opposite()->next()->next()->vertex()->id())
      beta = cot(pi-pn, pj-pn);
   else beta = 0.0;
   if (beta < epsilon) beta = 0.0;
   
   he->alphaplusbeta = alpha + beta;
}

void Mesh::floodfill_Dijkstra(unsigned idv) {
   
   cout<<"flood_fill(mesh) ... begin "<<endl;
   /// initializing 
   Vertex_iterator vi = this->vertices_begin();
   for ( unsigned int i=0; i<this->size_of_vertices(); i++) {
      vi[i].d = NINFINITY;
   }
   Vertex_handle h = vi[idv].halfedge()->vertex();
   
   h->d = 0;
   priority_queue<Vertex_handle, vector<Vertex_handle>, compareuv>   queue_vertices;
   queue_vertices.push(h);
   
   while (not queue_vertices.empty()) {
      Vertex_handle u = queue_vertices.top();
      queue_vertices.pop();
      Halfedge_handle hh = u->halfedge(), hhbegin = hh;
      do {
         Vertex_handle v = hh->opposite()->vertex();
         const Point3 &pv = v->point();
         const Point3 &pu = u->point();
         const Vector3 uv(pv[0]-pu[0],pv[1]-pu[1],pv[2]-pu[2]);
         const Number duv = sqrt(uv*uv);
         if (u->d + duv < v->d) {
            v->d = u->d + duv;
            v->previous(u);
            queue_vertices.push(v);
         }
         hh = hh->next()->opposite(); // 
      } while(hh != hhbegin);
   }
   cout<<"flood_fill(mesh) ... end"<<endl;
}

void Mesh::dijkstra_ShortestRing(const Vertex_handle &s, const Vertex_handle &t, list<Halfedge_handle> &ring) {
   /*
   priority_queue<Vertex_handle, vector<Vertex_handle>, compareuv>   Q, P;
   
   /// cleaning flags. It must be done before calling this function
   Vertex_iterator vi = this->vertices_begin();
   for (unsigned i=0; i<this->size_of_vertices(); i++) {
      vi[i].d = NINFINITY;
   }
   
   int ne=1; // number of vertices of the shortest path 
   const Point3 &pt = t->point();
   Vector3 m(pt[0],pt[1],pt[2]); // the centroid is initiallized with de source and target position
   
   /// The first step: the shortest path between s and t
   s->d = 0.0;
   Q.push(s);
   do {
      Vertex_handle u = Q.top();
      Q.pop();
      
      if (u == t) break;
      const Point3 &pu = u->point();
      Halfedge_handle hh = u->halfedge(), hhbegin = hh;
      do {
         Vertex_handle v = hh->opposite()->vertex();
         const Point3 &pv = v->point();
         const Vector3 uv(pv[0]-pu[0],pv[1]-pu[1],pv[2]-pu[2]);
         float duv = sqrt(uv*uv);
         if (u->d + duv < v->d) {
            v->d = u->d + duv;
            v->previous(u);
            v->previous_hh(hh);
            Q.push(v);
         }
         hh = hh->next()->opposite(); // 
      } while(hh != hhbegin);
   } while(not Q.empty());
   
   // setting all the vertices with d=oo
   vi = this->vertices_begin();
   for (unsigned i=0; i<this->size_of_vertices(); ++i) {
      vi[i].d = NINFINITY;
   }
   
   // Collect halfedges of first nearest path
   Vertex_handle tmp = t;
   while (tmp != s) {
      ring.push_front(tmp->previous_hh());
      const Point3 &pm = tmp->point();
      m = m + Vector3(pm[0],pm[1],pm[2]);
      tmp = tmp->previous();
      ++ne;
   }
   
   // Find the ring's center point
   {
   list<Halfedge_handle>::iterator itt = ring.begin();
   int cc = ring.size()/2;
   std::advance(itt, cc);
   const Point3 &pm = (*itt)->vertex()->point();
   m = Vector3(pm[0],pm[1],pm[2]);
   }
   
   typedef CGAL::Line_3<Kernel>                 Line;
   //~ typedef CGAL::Point_3<Kernel>                Point;
   //~ typedef CGAL::Vector_3<Kernel>               Vector;
   // computing the normal vector and the projected centroid
   //~ float *s_ = &this->vpositions[s->id()];
   //~ float *t_ = &this->vpositions[t->id()];
   //~ Point ss(s_[0],s_[1],s_[2]);
   //~ Point tt(t_[0],t_[1],t_[2]);
   
   const Point3 &ss = s->point();
   const Point3 &tt = t->point();
   Line l(ss,tt); // line betwwen the source and target
   Point3 pM = l.projection(Point3(m[0],m[1],m[2])); // punto medio projectado
   Point3 pm(pM[0],pM[1],pM[2]);
   //~ cout<<"pm projected on the line: "<<pm<<endl;
   //~ pm = ss + Vector3(pm[0],pm[1],pm[2]);
   const Vector3 n = Point3(m[0],m[1],m[2]) - pm; // normal vector
   Point3 pPlane = pm;
   Vector3 nPlane = n;
   
   //~ Vector3 on(n[1]-n[2],n[2]-n[0], n[0]-n[1]); /// the "on" Vector3 must be orthogonal to n
   //~ on = on/(sqrt(on*on)); /// The vector must be normalized because we need to compute the euclidean distance
   
   Vector3 on = cross_product(Vector3(Point3(tt[0],tt[1],tt[2]), Point3(ss[0],ss[1],ss[2])), n);
   on = on/(sqrt(on*on)); /// The vector must be normalized because we need to compute the euclidean distance
   
   //~ cout<<"ss: "<<ss<<endl;
   //~ cout<<"tt: "<<tt<<endl;
   //~ cout<<"pm: "<<pm<<endl;
   //~ cout<<"n: "<<n<<endl;
   
   float w = 1000.0;
   
   
   /// The second step: the opposite ring
   s->d = 0.0;
   P.push(s);
   do {
      Vertex_handle u = P.top();
      P.pop();
      
      //~ float *pu = &this->vpositions[u->id()];
      const Point3 &pu = u->point();
      if (u == t) break;
      if (u != s) {
         Vector3 pupm = Point3(pu[0],pu[1],pu[2]) - pm;
         if (n*pupm > 0) {
            u->d = 0.0;
            continue; // se estao do mesmo lado ---> cai fora
         }
      }
      
      Halfedge_handle hh = u->halfedge(), hhbegin = hh;
      do {
         Vertex_handle v = hh->opposite()->vertex();
         //~ float *pv = &this->vpositions[v->id()];
         const Point3 &pv = v->point();
         Vector3 pvpm = Point3(pv[0],pv[1],pv[2]) - pm;
         float don = on*pvpm;
         don = (don<0)? (-don):(don); // distance to plane defined by the normal vector "on" and the point "pm" 
         
         const Vector3 uv(pv[0]-pu[0],pv[1]-pu[1],pv[2]-pu[2]);
         float duv = sqrt(uv*uv) + w*sqrt(don*don);
         if (u->d + duv < v->d) {
            v->d = u->d + duv;
            v->previous(u);
            v->previous_hh(hh);
            P.push(v);
         }
         hh = hh->next()->opposite(); // 
      } while(hh != hhbegin);
   } while(not P.empty());
   
   // setting all the vertices with d=oo
   vi = this->vertices_begin();
   for (unsigned i=0; i<this->size_of_vertices(); i++) {
      vi[i].d = NINFINITY;
   }
   
   // Collect halfedges of second nearest path
   Vertex_handle tmp2 = t;
   while (tmp2 != s) {
      ring.push_back(tmp2->previous_hh());
      //~ float *pm = &this->vpositions[tmp2->id()];
      const Point3 &pm = tmp2->point();
      m = m + Vector3(pm[0],pm[1],pm[2]);
      tmp2 = tmp2->previous();
      ++ne;
   }
   */
}

void Mesh::save_closestBone(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ofstream fout(filename);
   if (fout) {
      fout<<"CLOSESTBONE"<<endl;
      fout<<_nvertices<<endl;
      Vertex_iterator vi = vertices_begin();
      while (vi != vertices_end()) {
         fout<<vi->rigging_group<<endl;
         ++vi;
      }
      fout.close();
   }
}

void Mesh::load_closestBone(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ifstream fin(filename);
   if (fin) {
      string info; uint nv; int id;
      fin>>info;
      fin>>nv;
      assert("CLOSESTBONE" == info);
      assert(nvertices() == nv);
      Vertex_iterator vi = vertices_begin();
      while (vi != vertices_end()) {
         fin>>id;
         assert(id != -1);
         vi->rigging_group = id;
         ++vi;
      }
      fin.close();
   }
}

void Mesh::save_closestJoint(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ofstream fout(filename);
   if (fout) {
      fout<<"CLOSESTJOINT"<<endl;
      fout<<_nvertices<<endl;
      Vertex_iterator vi = vertices_begin();
      while (vi != vertices_end()) {
         fout<<vi->rigging_group<<endl;
         ++vi;
      }
      fout.close();
   }
}

void Mesh::load_closestJoint(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ifstream fin(filename);
   if (fin) {
      string info; uint nv; int id;
      fin>>info;
      fin>>nv;
      assert("CLOSESTJOINT" == info);
      assert(nvertices() == nv);
      Vertex_iterator vi = vertices_begin();
      while (vi != vertices_end()) {
         fin>>id;
         assert(id != -1);
         vi->rigging_group = id;
         ++vi;
      }
      fin.close();
   }
}

void Mesh::save_skinningOnJoints(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ofstream fout(filename);
   if (fout) {
      fout<<"SKINNINGONJOINTS"<<endl;
      Vertex_iterator vi = vertices_begin();
      fout<<_nvertices<<" "<<vi->weights.size()<<endl;
      while (vi != vertices_end()) {
         for (uint i=0; i<vi->weights.size(); ++i) 
            fout<<vi->weights[i]<<" ";
         fout<<endl;
         ++vi;
      }
      fout.close();
   }
}

void Mesh::save_skinningOnJoints_assoc(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ofstream fout(filename);
   if (fout) {
      fout<<"ASSOC_JOINWEIGHT"<<endl;
      Vertex_iterator vi = vertices_begin();
      fout<<_nvertices<<" "<<0.0<<endl;
      while (vi != vertices_end()) {
         vector<Number> pesos; vector<int> ids;
         for (uint i=0; i<vi->weights.size(); ++i) {
            Number w = vi->weights[i];
            if (w > 0) { 
               pesos.push_back(w);
               ids.push_back(i);
            }
         }
         fout<<pesos.size()<<" ";
         for (uint i=0; i<pesos.size(); ++i) {
            fout<<ids[i]<<" "<<pesos[i]<<" ";
         }
         fout<<endl;
         ++vi;
      }
      fout.close();
   }
}

void Mesh::load_skinningOnJoints(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ifstream fin(filename);
   if (fin) {
      string info; uint nv; Number w; int njoints;
      fin>>info;
      fin>>nv>>njoints;
      assert("SKINNINGONJOINTS" == info);
      assert(nvertices() == nv);
      Vertex_iterator vi = vertices_begin();
      while (vi != vertices_end()) {
         vi->weights.resize(njoints, 0.0);
         for (uint i=0; i<njoints; ++i) {
            fin>>w;
            vi->weights[i] = w;
         }
         ++vi;
      }
      fin.close();
   }
}

void Mesh::save_skinningOnBones(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ofstream fout(filename);
   if (fout) {
      fout<<"SKINNINGONBONES"<<endl;
      Vertex_iterator vi = vertices_begin();
      fout<<_nvertices<<" "<<vi->weights.size()<<endl;
      while (vi != vertices_end()) {
         for (uint i=0; i<vi->weights.size(); ++i) 
            fout<<vi->weights[i]<<" ";
         fout<<endl;
         ++vi;
      }
      fout.close();
   }
}

void Mesh::load_skinningOnBones(const char *filename) {
   
   /// version using difussion over the mesh
   if (filename==0)                                return;
   if (_nvertices==0 or _nfaces==0)                return;
   
   ifstream fin(filename);
   if (fin) {
      string info; uint nv; Number w; int nbones;
      fin>>info;
      fin>>nv>>nbones;
      assert("SKINNINGONBONES" == info);
      assert(nvertices() == nv);
      Vertex_iterator vi = vertices_begin();
      while (vi != vertices_end()) {
         vi->weights.resize(nbones, 0.0);
         for (uint i=0; i<nbones; ++i) {
            fin>>w;
            vi->weights[i] = w;
         }
         ++vi;
      }
      fin.close();
   }
}
